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SUMMARY 


A viscous-inviscid interactive calculation procedure is developed for 
application to flow in cascades of two-dimensional airfoils. This 
procedure has essentially three components. First, a numerical solution of 
the Euler equations which can accommodate an arbitarily specified cascade 
geometry is carried out on a nonorthogonal curvilinear grid mesh that is 
fitted to the geometry of the cascade. A method of grid generation has 
been used which relies in part on a succession of conformal mappings. 
Second, a viscous solution for use in boundary layer and wake regions has 
been programmed. Finally, an interactive scheme which takes the form of a 
source-sink distribution along the blade surface and wake centerline is 
employed. Results have been obtained with this procedure for several 
cascade flow situations, and some comparisons with experiment are 


presented. 
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I, INTRODUCTION 

In recent years, a great deal of progress has been made in the 
development of faster, more efficient numerical procedures for the 
calculation of flow past aerodynamic shapes. Algorithms for the solution 
of the Euler equations and Navier-Stokes equations have been available for 
some time; for example, MacCormack's Method [1], an explicit time marching 
procedure, has been widely used since its introduction in 1969. More 
recently, time marching algorithms which have an implicit [2-4] or hybrid 
[5] character have been introduced. These implicit methods are not subject 
to the severe stability restrictions which explicit methods experience, and 
are therefore less time consuming. 

In addition to the need for an efficient flow calculation algorithm, 
another requirement in aerodynamic calculations is some technique for 
dealing with the complex geometries that often occur. Several different 
grid generation schemes have been developed in recent years to meet this 
requirement. Certainly among the most popular of these is a versatile 
method for dealing with aerodynamic geometries developed by Thompson, 
Thames, and Mas tin [6, 7]. In this method, which has been used in 
calculating the flow about isolated airfoils as well as other aerodynamic 
shapes, a non-orthogonal curvilinear grid mesh having a grid line 
coincident with the airfoil surface is generated by the solution of a 
system of elliptic partial differential equations. The coordinate 
transformation used in this method coresponds to the mapping of a region 
which encloses the airfoil in the physical plane, onto a region which is 
the interior of a rectangle in the computational plane. Steger [8] has 
combined the Beam and Warming implicit finite difference algorithm [4] with 
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Che grid generation procedure of Thompson et al, to simulate compressible 
flow about isolated airfoils. 

While similar in many respects to flow calculations for isolated 
airfoils, flow calculations for cascades encounter some additional 
difficulties in terras of the geometry and the boundary conditions which 
must be applied. The necessity for dealing effectively with complicated 
geometries in cascade flow problems has led to the development of several 
diverse geometry procedures (see for example [9—12 ] ) . Recently, Steger et 
al [13] have applied the approach used by Steger in the isolated airfoil 
problem to flow through cascades . 

In the present research effort, a body fitted nearly orthogonal 
curvilinear grid is generated by a method which relies in part on a 
succession of conform mappings. This is described in Section II. An 
implicit time marching finite difference solution of the Euler equations is 
then carried out on this grid in the manner described in References [8,13j, 
except for certain differences in the treatment of boundary conditions. 

The inviscid flow solution is discussed in Section III, In the present 
research effort, we have accounted for the effect of viscosity on the flow 
by coupling the inviscid calculation with a separate viscous shear layer 
calculation. This viscous calculation, which consists of a marching finite 
difference calculation for turbulent flow, is initiated at the stagnation 
point and proceeds through the blade boundary layers and into the blade 
wake. We have attempted to assess the economy of including viscous effects 
within the calculation in this manner, relative to the Navier-Stokes 
approach of Reference [13]. 
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II. GRID GENERATION 

The inviscid computations of the present work ere performed on a 
C-type body-fitted grid in which one family of lines forms open loops (C's) 
around the blade and wake. The grid is periodic and nearly orthogonal. 

This choice permits accurate resolution of the leading edge region and 
provides an appropriate location for the interactive wake boundary 
conditions to be described later. 

The grid generation employs two analytical mappings which take the 
multiply-connected exterior of a c^cade of airfoils to the interior of a 
simply connected domain. A numerically constructed mapping is then used to 
take this into a rectangular computational space. The first mapping 
transforms the exterior of a staggered cascade of semi-infinite flat plates 
in the z-plane into the interior of the unit circle in the w-plane. 

2 + A [a l°cjW - Ik) - ZC.0 sY |o^ c I'W)] (2-1) 

where is the stagger angle, s is the pitch, and (S/2.T0£X>p C LIT) 

This form is obtained from the standard mapping for a cascade of finite 
flat plates [14] by moving the singularities to 0 and +1 in the w-plane. 

At the leading edge z\ of the central plate dz/dw = 0. Solving for wj and 
substituting into Eq . (2-1) gives 

?o = 2,+^[VsmY + cosYlog(:UosY')] (2-2) 

When this mapping is applied to a real geometry, such as the turbine 
cascade in Figure 1, the flat plate is taken to run from a point just 
inside the leading edge through the downstream end of the wake. 

The second mapping transforms the interior of the unit crcle in the 
w-plane, with a branch cut from 0 to +1 along the real axis, to the 

interior of the infinite strip between the real axis and -iX/2 in the \ 

"3 -plane. 
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Note that reflection of ^ through the origin leaves w unchanged. This 
will be used later as a convenient means for analytical continuation. The 
image of the turbine cascade and straight wakes in the -plane (Figure 2) 
is a pair of parallel straight lines connected by a roughly s-shaped curve. 
Note the locations of corresponding points 1 to 5 along the contour in 
Figures 1 and 2 and the angle -fa between the flat plate and the wake. 
Reflection of "5, i.e. produces the opposite boundary of the 

analytically continued domain in the^-plane. 

In actual practice w is eliminated between Eqs . (2-1) and (2-3) and 
the transformation of the blade and wake from z to "3 is obtained by complex 
Newton iteration proceeding from point to point around the contour. To 
ensure that the branch cuts of the two logarithms are never crossed, the 
arguments of these logs are monitored and if either one changes by more 
than +71 between adjacent points, the value of the associated log at the 
new point is incremented by ^27Ii, i.e., in the opposite direction. 

The final mapping transforms the infinite strip in the "3-plane , 
bounded by the blade-wake contour and its reflection, into a rectangular 
domain with coordinates F = § + iV\_. If F is the complex potential for 
flow through the strip and we require F("3) = -F(-”5) and V\= -1 along the 


contour, then F can be represented as a contour integral. 

V OQ 


F(1) = ^ ^ W 1 °3 at + 


(2-4) 


- oO 


with C = (2/h)exp (-i/3). Here C, /3 , and h are, respectively, the complex 
velocity, flow angle, and channel width in the far field. Now set 
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I *1* t 


(2-5) 


where s is arc length and q^., q n are the vortex and source densities, 
respectively. We choose q n to cancel tv\e normal component of C. 



(2-6) 


This should minimize the magnitude of q£. For "5 on Lhe contour, IX* -1 and 


the imaginary part of eq . (2-4) gives an integral equation for q t . 


(2-7) 


xk S ll* Ha 0°g ?p|) 

= i + Ivn (Cl') 

Equation (2-7) is solved by a simple panel method with flat panels and 
locally constant q t and q n . Once q c is known, the 'ceal part of Eq . (2-4) 
gives 3 along the contour. 

oo 

VKe(cT)-^tV'Rfc(U \V£) f2 ' 8) 


oo 


%tTvn( cLs 


Generation of the grid in the rectangular <?,vo space proceeds in two 
stages. First, points are located on the boundaries such that the physical 
z-plane coordinates willbe periodic and continuous across the wake. As we 
shall see the ^\= 0 line transforms into the periodic boundary in the 
z-plane. Grid periodicity is enforced by distributing pairs of points 
symmetrically about the origin along the ^-axis . Continuity across the far 
wake is achieved by selecting a constant mesh spacing in this region such 
that the z-plane spacing is an integer fraction of (s)sin[/3 w | where $ is 
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the wake angle. Together these requirements imply that grid point location 
in is an odd function of § which becomes linear for large 11!. A 

fifth order polynomial is used for small l"$J in order to provide for 
clustering around the leading edge. The values of *5 along the contour are 

vs ^5 solution. In general the 
z-plane coordinates will still be discontinuous across the wake near the 
trailing edge as a result of contour curvature. A local straining is 
introduced to place a pair of points at the trailing edge and then pairs of 
neighboring points across the wake are adjusted until their z-plane images 
coincide. The distribution of points with V\at the two ends of the domain 
is arbitrary and a linear variation is used here. 

The interior values of are obtained from a finite difference 
solution of the complex Laplace equation 

q_ ^ i. - Q (2-9) 

* * v\> U 

Values of ^ are specified along V\= -1 and l^^^and the anti-symmetry 
property is used along 0. An ADI relaxation procedure is used to 

solve the finite difference equations with the Y\-inversion for fixed ^ 
performed simultaneously with chat for Estimates of the maximum and 

minimum eigenvalues of the ^ and Y\ matrices [15] are used to obtain a near 
optimum sequence of acceleration parameters. Figure 3 shows the grid 
distribution in the -plane for the turbine cascade of Figure 1. The 
upper plot boundary corresponds to V\» 0 and maps into the upper and lower 
periodic lines in the cascade plane. 






found by inverse interpolation in the *•$ 
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Figure 3. Turbine Cascade Grid 
in f-Plane. 
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Figure 4. Turbine Cascade Grid 


in z -Plane. 
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The final grid n the z-plane (Figure 4) is obtained by conformal 
mapping from the 5-plare using the two analytical functions (2-3) and (2-1). 
In this case continuity across the wake was obtained at the expense of a small 
amount of nonorthogonality. The rounded cap at the upstream boundary was 
obtained by extrapolation from the next two inner loops. Generation of this 
grid (99 x 7 points) required 1.4 seconds of CPU time on an IBM 3033 computer. 
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The Euler equations are written in a conservative form appropriate for 
a general curvilinear coordinate system. 


where 


Q = 


j. 

I 


IS. + jE + 3F_ Q 
at *3 + aw" ° 


(3-1) 


r * l 
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r /°u i 
/uU+V 
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■ ?V - 
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/° V 
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/OVU + ly-p- 

^ T 


e 




__ (e+f^V _ 


and where , 


T= 

V = V?y 

V= a n. x + v w,y 


In this equation, u and v are the x and y components of velocity, while U 
and V represent the contravariant velocity components in the ^ and V\. 
directions. These directions and velocity components are shown in Figure 5 
along with u s and u n , the physical velocity components that are 
respectively tangent and normal to an n,= const, grid line. Also,yO is the 
density, p is the pressure, e is the total energy per unit volume, and the 
fluid is assumed to obey the Perfect Gas Law. This equation was used by 
Steger in the study described in Reference [8], and since much of the 
theory of the i.nviscid solution used here is taken from that work, we will 
be content to present only a brief outline of this theory. 

The solution of equation (3-1) is accomplished by an implicit time 
marching algorithm, which is expressed in the "delta form" (see [4]) as, 
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Figure 5. The cascade coordinate system 
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- i^T^Xcr- o=-«(C*C) 
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N (3-2) 

/3g\ 

Here, A and B represent the Jacobian matrices Kq ) and l / > and the 
superscripts (n) and (n +1) indicate the time level at which a quantity is 
evaluated. and are the coefficients for the artificial dissipation 
terras , which have been added to the algorithm in both an implicit and 
explicit manner as suggested in [16]. 

Boundary Conditions 

For the description of the boundary conditions on the invisicid 
solution which follows, the reader may refer to Figure 5 where a typical 
computational grid mesh is displayed. Blade surface boundary conditions 
are obtained by first extrapolating ^ and U to the blade surface from the 
interior of the solution region. At the trailing edge, which is taken to 
be a cusp, extrapolated values of yO, p, and , are then averaged. The 

impermeability of the blade surface gives V = 0. Then the surface pressure 
is obtained by solving the tridiagonal system of equations that results 
from differencing the normal momentum equation 

= V vV-$ + 0< + vv*)fv 0-3) 

Boundary conditions along the wake centerline and periodic boundaries 
are obtained by averaging extrapolated values of /O, x>u, av, and e, at 
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coincident or periodic pairs of points. Along the unmatched portion of the 
wake centerline boundary, the averaging of extrapolated values is carried 
out between the unmatched boundary points and the point on the bottom of 
the centerline which is farthest downstream. This portion of the boundary 
of the computational domain is associated with the step-like form of the 
downstream boundary, which is introduced in an attempt to relieve the 
skewing in the coordinate system that would otherwise occyr in highly 
staggered cascades. It may be seen in Figure 5, and in the other 
computational grids displayed in this paper, that the coordinate skewing 
that results in these step-like grids is acceptably slight. The reduction 
in skewing is achieved at the expense of introducing this anomalous portion 
of the boundary; however, the treatment of this boundary segment did not 
prove to be a difficult problem since the calculation appeared to be 
relatively insensitve to the boundary conditions applied at this location. 

The treatment of the downstream boundary used here, follows the 
approach of Rudy and Strikwerda [17]. Their suggestion, for a rectangular 
computational domain with the downstream boundary oriented so that the 
outward normal is in the x-direction, was as follows. First, values of u, 
v, and T (temperature), are extrapolated to the boundary, and then the 
pressure is obtained by solving the equation, 

_/OC 'fe +0C (^~-M= 0 (3-4) 

where: 

c - local speed of sound 

p - exit pressure of the converged solution 
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OC - a parameter whose value Is chosen to optimize convergence 
By analogy, for the more complicated domain of Figure 5, we extrapolate 
values of u, v, and T, and obtain the pressure by solving, 





where , 


g' 


~Cix ) V const. ° r equivalently 



It may be seen in Figure 5 that the downstream boundary consists of two 
distinct and separate pieces; the approach to obtaining boundary values is 
the same at both the upper and lower portions of the downstream boundary. 

The solution of equation (3-4) is relatively easy to implement, since the 
only storage of values from the preceding time step required is at boundary 
points . 

At the upstream boundary, our approach is to specify v ; (?>) y and 
u + ^ c \. where X is the ratio of specific heats. A value for ( u - 2c \ 

' TTi' ^ tt' 

at each upstream boundary point is then obtained from 


Tt)-° (3- 

To solve equation (3-8), it is first rewritten as, 

(Tt +Cu '-°(^ IrV-o a- 

and then explicitly differenced in time. This simplified characteristics 
treatment of the upstream boundary is open to several objections, most 
notably perhaps in that the compatibility equation used here (equation 
(3-8)) is suitable for signal propagation in the negative x-direction. 



However, signal propagation in a direction that is essentially normal to 
the boundary is what on & - usually tries to accomodate with the method of 
characteristics treatment of boundaries (see for example (18]), and ic is 
easily seen in Figure 5 that the x-direction is at beat only somewhat 
normal to the upstream boundary. While a more sophisticated treatment of 
the upstream boundary would be more esthetically pleasing, we have chosen 
the present approach because it does not appear to degrade the numerical 
solution in any way and it is comparatively easy to implement. 

It is important to note that the boundary procedures described in this 
section were employed only after simpler approaches had failed. For 
example, the Rudy and Stridwerda approach to the downstream boundary was 
adopted after the method of extrapolating u, v, and T, and specifying p was 
found to cause the solution to become unstable. While such procedures are 
often employed in isolated airfoil calculations to assist with the 
convergence rate of the solution, it was our experience that in the absence 
of such precautions the cascade solution either converged to a result with 
noticeable errors or did not converge at all. It would appear that 
sophisticated boundary procedures, which can be an assistance for isolated 
airfoil calculations, are a necessity for cascades. In this regard we 
would mention specifically the nonreflecting downstream boundary condition, 
the characteristics treatment of the inflow boundary, and the use of the 
normal momentum equation to obtain the blade surface pressure. 
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XV, THE VISCOUS SOLUTION AND INTERACTIVE PROCEDURE 


The viscous solution consists of a finite difference marching 
calculation which is capable of dealing with blade boundary layers on both 
the pressure and suction surfaces, and of proceeding from there directly into 
the blade wake. The calculation can accomodate a flow which is compressible 
and turbulent; the turbulence modeling used is the two layer algebraic eddy 
viscosity of Cebeci and Smith [19], Additional features of the viscous 
solution include transition modeling [20], and a specialized turbulence 
modeling in the wake region [21]. The viscous calculation is initiated at 
the stagnation point in the following manner. First, the sign change of the 

velocity component u 3 is used to determine the two adjacent blade surface 

. . , . . . * 
points in the cascade grid which bracket the stagnation point . Second, 

the values of u s at these points are used to locate the stagnation point 

within the interval by interpolation, and the boundary layer is advanced to 

the bracket points by the similarity solution for stagnation point flow (see 

for example [22]). The suction surface and pressure surface calculations are 

then started from these locations and values. For the numerical marching 

calculation, we have employed a marching algorithm described in Reference 

[23]. 


Since the contravar iant velocity component U might seem to be the 
most obvious quantity to monitor for a sign change to locate the stagnation 
point, it is noted in passing that u s and U do not necessarily share the 
same sign at blade surface points if there is a non-zero velocity component 
normal to the surface. A surface injection distribution is specified for 
subsequent inviscid solutions within the viscous-inviscid iterative 
process, so a non-zero normal velocity component can be expected in this 
calcuation. In this case, u s and not U is the appropriate choice. 
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Governing Equations 

For the viscous marching calculation we take as the governing equations; 


'L.(/ oLL ') 4 L(/> v ) - n 


(4-1) 


/OUL^iis 


+ + -/(mV)) 


V ' ^ y 


(4-2) 



(4-3) 


+ 1 _ (AL ih 

V vPv 1y 


~/o(hV)) + (^ 



"Z 2 — /° ft T 


(4-4) 


V> = c P T 

where Pr, Cp,andto, are taken as constants, and the subscript oo denotes 
evaluation at the free stream. (In contrast to the notation used with the 
inviscid solution, x and y are now used to denote the streamwise and 
streamnormal directions.) Furthermore if we write, 


(4-5) 

(4-6) 


-/° ( u-' V -) = a t 



(4-7) 
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and , 




and let 


4 + ^ 

Cm* "pT j 


4 + = lii 


V 


(4-8) 


, and Pr t ■ , where Pr^. is a 

specified constant, then Equations (4-2) and (4-3) may be rewritten as; 


/ ou -'^ + 


(/° V % =/*•“*>% + \y( M ^ 0 \) <*-»> 




+ V^(' + Cfe^v) (4 ‘ io) 


Equations (4-1) and (4-4) - (4-10) form the basic system of equations which 
the marching calculation solves, exclusive of the turbulence modeling. To 
model the turbulence, the two-layer algebraic eddy viscosity of Cebeci and 
Smith [19] has been used. 

For y<yi , £ m =k m [ 
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L- o.s y(l- exp(-y//C)) 

K = itO’/u^V^’/N 


(4-11) 


(4-12) 


(4-13) 
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(4-16) 


(4-17) 


where the subscript w indicates that a quantity is evaluated at the wall 
(y » 0 ). 

R~ r y > yi » £m ~ ^ m 0 


^ = 0.0168 
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(4-18) 


The boundary between the inner and outer layers (y = y^) is taken as the 
location where • 

Equations (4-11) - (4-18) are used to treat the blade boundary layers. 
In the wake the modeling used is somewhat different, and follows the 
three-layer approach described in Reference [21]. 

For y<y 2 , ~ ^ 


O.y (^x) TE Y { 


(4-19) 
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where y “ y c when u a 1.10 (u)x,E.> and the subscript T.E. denotes evaluation 
at the trailing edge. 

g°r„,y.2..^y...<yi > £ . TO .* fen - t 

Equation (4-11) is applied in this region with, 


L Y (4-20) 

ESLJ.> yj A. m a Axt i 0 

Equation (4-18) is applied in this region without alteration. 

The boundaries between the layers, y » yj_ and y “ y£> are taken as the 

locations where and £va"£wV ’ res P ect *- ve ly • Furthermore, if 

( - . \ A , then t^is taken to equal A for all y. 

V yv v s ^ '-vvio * o 

An additional feature of the turbulence modeling used in this study is the 
transition intermi ttency factor (X tr ). This Y tr , which multiplies the eddy 
viscosity determined from the preceding formulas, is calculated according to 
Equations (4-21) and (4-22) which follow from the work of Chen and Thyson 
[2GJ. 


l-exp(-G( (j ^oLx) 
& - .000 


l? *■ 
co 


ft 


1.34 


(4-21) 


(4-22) 


A value for the transition Reynolds number (Re x ) is specified for the 

tv- 


calculation. 
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Grid Meah for the Marching Calculation 

The viscous inarching calculation is carried out on a grid that io 
considerably more dense than the one used for the inviscid solution. In the 
streamwise direction along the blade surface and wake centerline, 30 boundary 
layer grid divisions correspond to each single inviscid solution grid space. 

In the direction normal to the blade and wake, a variable grid spacing 
technique described in Reference [19] has been employed in an attempt to 
better resolve the large gradients which characterize a turbulent shear layer. 
For this approach, in wnich the ratio (K.) of consecutive grid spacings is 

a K) , values for the first grid interval (Ay)j[ and the 

j 

ratio (K) are specified. 


fixed (i.e., 


Ca/) 


The Numerical Scheme 

As mentioned previously, a marching algorithm described in [23] has been 
used in the viscous solution. For the grid system described in the preceding 
section, Equations (4-1), (4-9), and (4-10) result in the following finite 
difference equations. 


/Pj+ijj+i _ (/ Q V)[f..^r(/ , v) 

o>o u. 


(4-23) 


a, u; . u r,U~, LL ^i . u , 

d (AKV W , +( ^ v, m c / Vl u vi UxVi+1 

+ ^k y jCIK«Vn)i,j)^ 

V (AyV; O + Kl A 1 + '■'*'> ) \ CAy)j O+K) / 

... . /(^mV/i+rC^ljj-A/ ut+vuru-'m, j-i N 

M bi aryijW-Jv TtyYiOtW'J 


(4-24) 


24 


offlumi pi\zz its 

OF POOR QUALITY 


/i. .ji, . V\ t+ 1 ^ \ ~Vu)j /*-Ta 

AdN + y v > i j (W+k') 

i . . . / ( „ \ (h^L+r(Vi«>V\N / 1 f \ \ AM* 


VMjO 

+ iiij f . + //+\ A Vhui,.urOuq\v l>ll ] A*t„j-i\ -L 

+ tv ^H^Pv-A (KCK+ft/aX*^ / 

I / A WI A ^j-lV I i/i A VwwrfN x 

Pr- t C^jd+KVV 1 Rr t a Ay^ (H-KV 


(Aft. (I+W) 


(4-25) 


(Ay)jO + K) 


The system of difference equations consisting of (4-23), (4-24), and 
(4-25), together with the various property relations and turbulence 
equations, are solved subject to appropriate boundary conditions These 
boundary conditions are applied at free stream, blade surface, and wake 
centerline locations; it is in this regard that we mention the one-sided 
difference representation of cr. 


^ Vi _ - V\uj»a.+ H Vn,jM~3 ^ Ai 
V L,j“ 0-K)(uy) J+l 


(4-26) 


The procedure for advancing the solution from the (i) station (presumed known 
to the (i+l) station (presumed unknown) is as follows: 

(i) Solve Equation (4-24) to obtain u^+ij for all j. A tridiagonal 
inversion is required, 

(ii) Solve Equation (4-25) to obtain h^ + j_ j for all j. A tridiagonal 

inversion is required. (Due to the adiabatic wall and wake center- 

\L 

line symmetry boundary conditions, ® at y = 0, the matrix is 

not tridiagonal, but it will take this form after one Gaussian 
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elimination step. A similar preconditioning of the matrix is 
required when solving Equation (4-24) in the wake). 

(iii) Equation (4~4) and (4-5) and (4-6) are solved for t j , y^i+1 } j > 

and/li+ij. 

(iv) Equation (4-23) is marched from y 3 0 to the outer edge of the 
shear layer to obtain (yOv )•£+]. . 

(v) The turbulence modeling equations ((A-ll)- (4~22)) are solved to 
obtain (£ m + )i+i f j. 

The solution has now been advanced to station (i+1). Although this numerical 
scheme is a low order method, we have found it to be simple, inexpensive, and 
very stable. We have perceived these virtues of the method in a previous 
study [24], as well as in the present work. The reliability of the method 
was an important feature in a calculation that was expected to encounter 
laminar-turbulent transition, boundary layer separation, and the sudden 
change in boundary conditions that occurs when marching off the trailing edge 
and into the wake. It was for this reason that our first choice of a 
inarching algorithm, the Keller Box Method (see for example [19]), which is 
significantly more accurate than the present scheme, was replaced with a 
method which appeared to be less sensitive. We have attempted to compensate 
for the lower accuracy of the method by an increase in grid density. 
Additional Computational Features 

Two additional computational features of the viscous solution remain to 
be discussed. As the first of these, it is noted that the calculation 
procedure includes provisions for extending the grid in the y-direction, 
should the shear layer approach the outer edge of the grid too closely. This 
feature enables the calculation region to grow with the boundary layer and to 
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contain it as economically as possible. Second, to equip this direct 
boundary layer solver with at least some ability to cope with regions of 
separated flow, we have altered the streamwise convection terms in Equations 
(4-24) and (4-25) for reverse flow in the manner of Rehyner and Flugge-Lot 7 ; 
[25]. 

The Interactive Procedure 

The interactive procedure in this calculation consists of an iteration 
between the inviscid and viscous solutions previously described. The effect 
of the presence of the viscous shear layer on the inviscid solution, is 
modeled as a source-sink distribution along the blade surface and wake 
centerline. This source-sink distribution is obtained from the viscous 
solution according to the following expression; 



co 



(4-27) 


where the notation used in the viscous solution (Section IV) appears on the 
right side of this equation, and the left side conforms to the notation of 
the inviscid solution (Section III). Once values of ^u n ) have been 
calculated and supplied to the succeeding inviscid solution, there are some 
alterations in the treatment of the blade surface and wake centerline 
boundary conditions which we described Section III, 

Blade surface boundary conditions retain the form described earlier, 
except that the specification of V = 0 is replaced with the specification of 
(yOu n ). The surface pressure is still obtained by a tridiagonal numerical 
solution of Equation (3-3). Along the wake centerline, values of (yOu n ) 
calculated from the preceding viscous solution are summed at coincident 


points, and these sums are taken to represent a discontinuity in the values 
of }y3u n [ at those locations. The equality of y3,u s , and p, at coincident 
points is enforced by averaging extrapolates. 
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V. NUMERICAL RESULTS AND DISCUSSION 


The results that have been obtained with the numerical procedure 
described in the preceding sections are presented here. The procedure has 
been applied to three cascades: an unstaggered cascade of NACA 0012 blades 

with a gap-chord ratio of (1.0), a 45° staggered cascade of NACA 65-410 
blades with a gap-chord ratio of (0.777860), and a 28.5° staggered cascade 
of NACA 65 — (12 ) 10 blades with a gap-chord ratio of (1.0). Results are 
presented for several different flow situations in each of the three 
cascades including different Mach numbers and angles of attack. In certain 
cases where the viscous inarching solution was incapable of dealing with the 
separated regions that occurred, and it was therefore impossible to perform 
an interactive calculation, we have presented results from the sole 
inviscid solution. The cascade of NACA 0012 blades and the cascade of NACA 
65-410 blades were chosen as simple test cases to evaluate the performance 
of the numerical procedure. The cascade of NACA 65-(12)10 blades has been 
the subject of an experimental study of Briggs [26], and was chosen as it 
afforded the oportunity for comparison with experimental results. It is 
hereafter referred to as the Briggs cascade. The results obtained in each 
of the three cascades are discussed separately. However, the blade chord 
length is .25 ft and the flow is air in all of the test cases considered. 
Unstaggered Cascade (NACA 0012 blades) 

An unstaggered cascade of NACA 0012 blades with a gap-chord ratio of 
(1.0) was chosen for use in the initial tests of the numerical procedure, 
with the idea that an unstaggered cascade of symmetric airfoils would 
minimize the geometry related difficulties and would permit attention to be 
focused on other computational aspects of the method. The grid system used 
for this cascade is shown in Figure 6. 
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Figure 6. Grid system, unstaggered cascade 
( N A C A 0012 blades). 
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A subsonic flow at zero angle of attack was chosen as a first test 
case for the interactive calculation procedure. More precisely, for the 
upstream and downstream boundary conditions on the invscid solution we have 
specified the following: 


Upstream Tangential Velocity (vj ) - 0. ft/sec 

X~ — °o 


Upstream Riemann Invariant C( u+ ~ 

5980. ft/sec 

Upstream Isentropic Constant 

948210, (p - Ibf/ft^yO- slugs/ft^) 

Exit Pressure of the Converged Solution “Pfeo “ 2125. lbf/ft^ 

The specifications result in a velocity of 496 ft/sec and a Mach number of 
(.45) at the upstream boundary. 

The interactive calculation procedure was run for four global 
viscous-inviscid iterations in this test case with no apparent difficulty; 


a sampling of the results may be found in Figures 7-9. In Figure 7, values 

- •p~'p° o 


of the pressure coefficient on the blade surface 


(C p = 


from the first and last (fourth) inviscid solutions are displayed. Surface 


S- 00 ^ 


pressure coefficient values have been plotted for all grid points on the 
top and bottom surfaces of the blade; however, due to the symmetry in this 
flow situation, these values for corresponding points on the top and bottom 
are largely indistinguishable. This figure demonstrates the expected 
result that for the case under consideration, the fluid viscosity has only 
a slight influence on the surface pressure. In Figure 8, the convergence 
history of each of the four inviscid solutions is shown. The maximum 
residual, the decay of which is monitored in these plots, is calculated 


according to the expression ^ ^ 


oo 


f\r\* 


where &p is the change in 
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Figure 7. Surface pressure coefficient 

(NACA 0012 cascade, M!=. 45) 









Figure 9. Mass injection rate* 


pressure that occurs between two consecutive time steps, p^ is the exit 
presssure, and where points along the blade surface and wake centerline 
have been examined. The trend towards convergence, in both the time 
marching inviscid solutions and the global viscous-inviscid iteration, is 
apparent. Further evidence of the convergence of the method is shown in 
Figure 9, where plots showing the behavior of the mass injection rate (^u n ) 
for each of: the four global iterations are found. This behavior is 
examined in a region centered about the trailing edge (t.e.), as this 
region is one of special importance and sensitivity. Again, a trend 
towards convergence is indicated. 

For this calculation, a time step was used which was estimated to be 
about 20 times larger than the Courant- Friederichs - Lewy (C.F.L.) Limit 
based on the smallest grid spacing in the field. This value was chosen as 
it was about the largest time step that could be used without destabilizing 
the calculation. The time step was held fixed throughout the calculation; 
no attempt was made to improve the convergence rate by cycling or varying 
the time step in any way, nor was it necessary to use a smaller initial 
time step to accomodate the calculation's impulsive start. (The only 
concession made to this impulsive start was to enforce the impermeability 
of the blade surface gradually, over 50 time steps.) The time step size 
did not differ greatly from this value in any of the test cases described 
in this section. For the artificial dissipation terms that appear in the 
inviscid solution algorithm, values of 6^- *05 and (£j_/(: e ) = 2 
were used in this test case, and similar values were used in all subsequent 
calculations. The value of 2 followed from suggestions made in Reference 
[16], while £r e = .05 was chosen as it was about the largest value which did 
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not destabilize the calculation. The numerical evidence of this and 
subsequent calculations was that a stability limit of approximately (.06) 
existed for (r e . This observation conflicts with Reference [16], where the 
introduction of implicit damping (£i) is advocated as a means of relieving 

a stability limit on the basic algorithm (i.e. = 0) of (** 0.625). 

I(r> 

This point will receive additional comment later in this section. 

As a second test case for the the interactive calculation procedure, a 
flow into the unstaggered cascade at a 2° angle of attack was considered, 
with a Mach number of (.45) at the upstream boundary. This flow situation 
was chosen as a simple test of the ability of the method to treat a lifting 
cascade. An interactive solution, of three viscous-inviscid iterations was 
performed, the convergence history of which differed very little from the 
preceding test case. Values of the surface pressure coefficient from the 
first and last inviscii solutions are shown in Figure 10. Again, the fluid 
viscosity has an effect on the surface pressure distribution which though 
noticeable is small. 

In the third and fourth test cases, situations have been considered 
where the effect of viscosity on the flow was more pronounced and the need 
for a viscous-inviscid interactive capability more clearly demonstrated. 

Both of these test cases involve flows in the unstaggered cascade at zero 
angle of attack, but higher Mach numbers. The third test case involved a 
subsonic flow with an upstream Mach number of (.63). The surface pressure 
coefficient for this case, which is shown in Figure 11, displays a larger 
difference between the first and last solutions in the iteration, than 
appeared in the previous test cases. The effect of the viscosity on the 
surface pressure is still more apparent if we consider a transonic flow 
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Figure 10. Surface pressure coefficient 

(NACA 0012 cascade, M|~4 
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Figure 11. Surface pressure coefficient 

(NACA 0012 cascade, M|=. 
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Figure 12. 
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Figure 13. Maximum residual history. 
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situation such as the fourth test case, with an upstream Mach number of 
(.65). The surface pressure coefficient for this test case, shown in 
Figure 12, indicates that a shock is now present in the cascade. While the 
shock resolution in this case appears to be only fair, this shock smearing 
enables our direct boundary layer solver to negotiate the region of shock 
irnpingment on the boundary layer, and makes a viscous-inviscid iterative 
solution possible in this test case. We will return to the discussion of 
shock resolution and its implications for the viscous solution later in 
this section. For the present it will suffice to mention than the 
iterative solution converged successfully although not as quickly as in the 
subsonic test cases, and that the convergence history of this solution is 
found in Figure 13. 

45° Staggered Cascade (NACA 65-410 blades) 

A cascade of NACA 65-410 blades staggered at a 45° angle with a 
gap-chord ratio of (0.777860) was chosen to test the numerical procedure on 
a more geometrically difficult cascade. The grid system used for this 
cascade is shown in Figure 14. A subsonic flow that was well aligned with 
the blades was chosen for the first test case, with a Mach number of (.47), 
and a flow angle of 48.4° (i.e., an angle attack of 3.4°) at the upstream 
boundary. An interactive solution of three iterations was performed, the 
results of which are displayed in Figure 15. Also, the convergence history 
of this solution is shown in Figure 16. 

An interesting feature of the surface pressure distribution in Figure 
15 is the very low pressure that occurs on the pressure surface near the 
leading edge. In the second test case for this cascade, we sought to 
eliminate this feature by increasing the flow angle at the entrance. The 
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Figure 14. Grid system, 45°staggered cascade 

(NACA 65-410 blades). 
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Figure 15. Surface pressure coefficient 

(NACA 65-410 cascade, M = 
3.4' angle of attack). X 
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Figure 17. Surface pressure coefficient 

(NACA 65-410 cascade, M|=79, 
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Figure 18. Surface pressure coefficient 

(NACA 65-410 cascade, Mj = 

'X 

7.3'angle of attack). 




speed of the flow was also increased so that the flow at the entrance 
exhibited a Mach number of (.79), and a flow angle of 52.5° (7.5° angle of 
attack). It may be seen in Figure 17 that the increase in angle of attack 
has effectively removed the low pressure region on the pressure surface. 
Also, the figure indicates the presence of an incipient shock on the 
suction surface at about 25% of chord, due to the increased speed of the 
flow. This shock is more apparent if the speed of the flow is increased 
slightly. This has been done in the third test case, where the entering 
flow had a Mach number of (.83), and a flow angle of 52.3°. As the surface 
pressure coefficient distribution shown in Figure 18 indicates, a shock has 
now formed, the location of which is clearly dependent on viscous effects, 
While convergence tended to be somewhat slower for test cases with the 
staggered, cambered blade cascade than those with the unstaggered cascade 
of symmetric blades, nevertheless the interactive calculation procedure was 
judged to have performed well in all cases considered. 

28.5° Staggered Cascade (NACA 65-(12)10 blades) 

The success that was experienced in a variety of preliminary Lest 
cases prompted the application of the numerical procedure to a cascade 
situation where a comparison with experiment could be made. The cascade 
chosen for this purpose consisted of NACA 65-(12)10 blades at a stagger 
angle of 28.5° and spaced with a gap-chord ratio of (1.0). The grid system 
used for calculations in this cascade was displayed previously, in a 
different context, in Figure 5. Comparisons have been made with data 
collected by Briggs [26] over a range of Mach numbers. In these 
comparisons we have restricted our attention to the data in [26] which 
satisfies that author's stated two-dimensionality criterion. The stated 
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flow angle in all Sriggr' data is 45°, which corresponds to an angle of 
attack of 16.5°. 

A somewhat disappointing result of the numerical calculations carried 
out on the Briggs cascade, was that a significant boundary layer separation 
occurred on the suction surface near the trailing edge, which caused the 
viscous inarching solution to fail. This boundary layer separation 
prevented the execution of viscous-inviscid iterative solutions; 
consequently, the numerical results presented for the Briggs cascade are 
strictly inviscid resuts . The inviscid procedure did however perform quite 
well in these calculations. 

We first consider two subsonic test cases with upstream Mach numbers 
of (.42) and (.61). 


/ (s _ I X — - on \ 

Values of the surface pressure coefficient l \j D — ”1 _ , , a» ) 

\ P al/ 3 U | x --oo / 

from the numerical solution and experiment are compared for the .42 

case, in Figure 19; the agreement apears to be good. Also, the pressure 


rise predicted by the numerical solution, 


1.055, agrees well 


"P | x =>-<*> 

with the experimental value of 1.055. For the M | x - = .61 case, a 

comparison of surface pressure values in Figure 20 shows fairly good 
agreement between the numerical solution and experiment. Also, the 
calculated value for the pressure rise in the cascade of 1.148 agrees 
fairly well with the experimental value of 1.135. 

While the agreement between these two pressure rise values is 
acceptable, it is not as good as in the M | X=S _ O 0 = .42 case, and the 
comparison deteriorates somewhat further in the transonic case that follows. 
In that transonic test case, for which additional results will be presented 
shortly, the apparent overprediction of the pressure rise is still more 
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Figure 20. Surface pressure coefficient 
(Briggs cascade, M|=.61). 
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Figure 21. Surface pressure coefficient 
(Briggs cascade, MI—61, 

. lx = -oo 

42.8 flow angle). 







1 


severe. It is reasonable to ascribe some of the disparity between the 
numerical and experimental results to the absence within these numerical 
solution of any provision for viscous effects. This reasoning is supported 
both by the fact that the numerical procedure overpredicts the pressure 
rise, and by the observation that this overprediction tends to worsen with 
the increase in Mach number. However, our experience with interactive 
solutions carried out in preliminary test cases indicated that it was 
unreasonable to attribute all of the overprediction to viscous effects. 

For this reason, it was resolved to investigate the possibility that the 
effective angle of attack in the experiment in Reference [26] was slightly 
less than the stated value. It was also judged a possibliity that 
the pressure rises recorded in [26] were more reliable data than the flow 
angle. Our approach then to the investigation of this possibiity was to 
recalculate the second test case with the 45° flow angle requirement 
removed and a smaller value of the pressure rise enforced. This 
recalculation resulted in a Mach number of (.61), but a pressure rise of 
1.138 which was closer to the experimental value. The numerical solution 
gave a value for the upstream flow angle of 42.8° (i.e. a 14.3° angle 
of attack). The surface pressure coefficient plotted in Figure 21 appears 
to be in slightly better agreement with the experimental data than the 
previous numerical solution in Figure 20. 

As the final test case wth the Briggs cascade, we have considered a 
transonic flow situation at an upstream Mach number of (.76), and a flow 
angle of 45°. 

The numerically generated surface pressure coefficient is compared 
with the expe7.‘iraental data in Figure 22. A value for the cascade pressure 
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rise of 1.265 was calculated, the measured value being 1.218. Aside from 
the shock oscillations in Figure 22, it is clear that the agreement between 
the numerical and experimental results is not particularly good. The 
discussion of he shock oscillations is postponed until later. The general 
lack of agreement in Figure 22 together with the disparity in pressure rise 
values however, suggested that recalculation of this test case based on a 
matched pressure increase rather than a matched flow angle, was in order. 
Since it was unclear as to what portion of the disparity in the pressure 
increase should be attributed to viscous effects, we have performed two 
such recalculations at different pressure rise values In the first of 
these cases a flow with m| x = -oq" *76, and a pressure rise of 1.248, 
resulted in a flow angle of 42.5° (i.e. a 14° angle of attack). In the 
second a flow with M| x = - 03 = .76, and a pressure rise of 1.217, resulted in 
a flow angle of 39.2° (i.e. an angle of attack of 10.7°). The numerical 
results are cempared with the experimental data in Figures 23 and 24. 

Of the comparisons in Figures 22-24, perhaps the best agreement is in 
Figure 23, although the agreement in none of these solutions is especially 
good and the appearance of disagreement is further increased by the 
spurious shock related wiggles that are present in these plots. In closing 
the discussion of the present test case, it is noted that the numerically 
generated pressure coefficient distributions in Figuns 22-24 demonstrate 
that a relatively small change in the angle of attack can cause a large 
change in results. Since differences in the effective angle of attack of 
this magnitude might well fall within the bounds of experimental error, the 
comparisons using matched pressure rises are to be preferred. 
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Discussion 


The numerical results presented here demonstrate the applicability of 
this numerical calculation procedure to the analysis of flow in cascades. 
The calculations for which results have been presented required from one to 
three hours (C.P.U. time) on a CDC Cyber 170/760, where most of the expense 
was associated with the inviscid solution. Although these run times are 
large, they were not viewed as prohibitive, particularly since faster 
computing machines are available. While the method in, its present form has 
been successfully applied to a number of cascade flows, our experience with 
the method indicates that its range of application could be greatly 
extended with some modifications. It is our perception that the two most 
important limitations on the method at present are the shock resolution 
problems in the inviscid solution, and the failure of the viscous marching 
solution (and consequently the viscous-inviscid iteration) in situations 
with significant separation. 

The presence of shock induced oscillations in the transonic test cases 
constitutes the most serious shock related difficulty encountered with the 
method. These oscillations are prominently displayed in Figure 22, for 
example. Some standard remedies, which have been applied to this problem 
in isolated airfoil calculations, proved incapable of relieving the problem 
for the cascade flows considered here. For example, the attempt to 
suppress these oscillations by Increasing the damping coefficients was 
frustrated by an apparent stability bound for this calculation of 
Also, an attempt to remove the oscillations by the introduction of 
"conservative spatial switching", described in Reference [4], and by 
transitioning the switching operator as described in [4], failed when the 
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calculation became unstable in the presence of supersonic flow. The 
failure was somewhat surprising in view of the success experienced with 
this approach in previous studies [4,8], In our attempt to resolve the 
problem in thi" manner, we have dropped the damping terms in supersonic 
regions and have in general tried to conform as closely as possible to the 
prescriptions of [4,8], with the sole exception that we employed a lower 
order upwind differencing. 

Both the failure of the attempt to incorporate "conservative spatial 
switching" wiohin the method, and the apparent stability bound of 
seemed to indicated that flow in cascades was in some way a more severe 
test of the inviscid solution algorithm than flow past an isolated airfoil. 
It is possible that these failures occured because of some incompatibility 
of the larger damping and the switched differencing with boundary 
procedures used in the inviscid solution. Also, stability analyses of this 
algorithm (for example [16]) suggest that more restrictive stability bounds 
result from the application of periodic boundary conditions, but these 
analyses do not predict bounds as restrictive as our numerical experience 
with the present method would indicate. Finally, the possiblity always 
exists that an error in the coding remained undetected, although it is 
unlikely that this was the source of the problem. Regardless of the 
reasons for the failure of these attempts at smoothing the shock 
oscillations, it is clear that this difficulty must be overcome for the 
method to be applied with confidence to shocked flows. 

The second important limitation on the method was the failure of the 
viscous calculation in situations with significant separation. Unlike the 
difficulty in removing the shock oscillations, this problem was not 


CjL 


unanticipated. The present finite-difference viscous marching solution was 
chosen for its simplicity, accuracy, and reliability in the treatment of 
attached boundary layers and wakes. However, it was not expected that this 
march!. g solution would accomodate flows that were severely separated. 

While the viscous marching solution currently employed performed well in a 
variety of cascade flow situations, the disappointing failure in the Briggs 
cascade indicated the desirablity of some improvement with regard to 
separation. As this would likely require the use of an inverse boundary 
layer procedure, which would then be coupled with the present direct 
inviscid solution, a viscous-inviscid iterative scheme of the type referred 
to as "semi-inverse" (see for example [27]) would seem to be required. 

In conclusion, it is our opinion that the interactive calculation 
procedure developed in this study constitutes a useful tool for the 
analysis of cascade flows which are unshocked and only mildly separated; 
and with some modification, the generality of the procedure could be 
increased with regard to these features. The procedure is quite general in 
terms of cascade geometry and can accomodate a wide range of blade shapes , 
blade stagger angles, and blade s pacings. While a Navier-Stokes approach 
to the inclusion of viscous effects within the numerical calculation would 
undoubtedly be required in certain severe flow situations (e.g. massive 
boundary layer separation with vortex shedding), the viscous-inviscid 
interactive approach provides an alternative in the analysis of less severe 
flow situations, such as a cascade operating at or near design conditions. 
In a time marching solution of the Navier-Stokes equations, it is unlikely 
that the computational effort expended per grid point would greatly exceed 
the value for a time marching solution of the Euler equations. However, it 
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is an advantage of the interactive approach, that this inviscid time 
marching solution is carried out on a grid that is sparse in comparison 
with the grid requirements of a Navier-Stokes solution. 
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